Volumetric formulation of lattice Boltzmann models with energy conservation 
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We analyze a volumetric formulation of lattice Boltzmann for compressible thermal fluid flows. 
The velocity set is chosen with the desired accuracy, based on the Gauss-Hermite quadrature proce- 
dure, and tested against controlled problems in bounded and unbounded fluids. The method allows 
the simulation of thermohydrodyamical problems without the need to preserve the exact space- filling 
nature of the velocity set, but still ensuring the exact conservation laws for density, momentum and 
| energy. Issues related to boundary condition problems and improvements based on grid refinement 

^v^j . are also investigated. 
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Recent studies on the lattice Boltzmann method (LBM) p], Q have prompted tremendous advancements in the 
capabilities of the method to systematically handle and reproduce complex flow properties (3l4l3j. When dealing 
with isothermal Navier- Stokes equations with small degree of compressibility, LBM is frequently used with standard 
lattices possessing a relatively small number of velocities (less than ten in two dimensions and just a few tens in 
three dimensions). However, the situation is quite different and deeply challenging for compressible thermal flows 
EMEUS- As a matter of fact, it is not easy to incorporate the temperature into the lattice equilibrium when using 
the standard lattices, and simultaneously to satisfy a number of conditions for recovering the correct thermo hydro- 
dynamical description of compressible flows. This has triggered the development of higher order LBM schemes with 
larger and more isotropic sets of velocities [ll|, [TH, [H| • Possible ways of obtaining these models are by discretizing the 
Boltzmann equation on the roots of Hermite polynomials and systematically derive new complete Galilean-invariant 
LBM schemes [ll|, [HI, [l8, 20-23], or also introduce a systematic approach to construct higher-order lattices for stable 
LBM based on the entropic approach [l2|, [24| . 

Whatever is the systematic procedure used, when the roots of the velocities are irrational, the corresponding discrete 
velocities cannot be fitted into a regular space-filling lattice. Thus, one of the most important advantages of the 
LBM, i.e. the exact space discretization of the advection step, is lost for the off-lattice models. For achieving a better 
accuracy, still keeping exact space-filling discretization, LBM with a large number of velocities were suggested based 
on the Hermite-Gauss quadrature procedure. Just to give an example, the D2Q53 and D2QS1 models detailed in a 
recent number of papers [H, [H, [25| , allow for a precise higher order accuracy, but they possess much less flexibility 
with respect to standard models due to the increasing number of kinetic fields. Also, the use of those models with an 
exact discretization of the streaming step may pose the serious problem of boundary conditions, which is not an easy 
task when the number of velocities is increasing. 

In order to keep a reasonably high accuracy of the lattice velocities and still retain a tractable number of them, one 
is somehow forced to move on off-grid lattices and find the correct computational scheme to be used. In particular, 
these have included interpolation schemes [26[ , different finite volume schemes [l|, [27H29] and LBM with local grid 
refinements and unstructured grids/adaptive meshes 1 30, 31]. A particularly interesting approach has been discussed 
in a recent number of papers by Peng and coworkers [27H29|, based on finite volume techniques in the LBM framework. 
The resulting lattice Boltzmann schemes integrate the differential form of LBM using a finite- volume scheme in which 
the unknown populations are placed at the nodes of the mesh and evolve based on the fluxes crossing the edges of the 
corresponding elements. 

In this paper we numerically and theoretically explore the potentiality of a volumetric formulation for LBM with 
active thermal fluctuations (hereafter refereed as TVLBM). The thermal part of the model heavily relies on Hermite 
quadratures with integer and non integer roots [ll|, [H, [32| for both bounded and unbounded flows. This kind of 
approach has the obvious disadvantage to loose the exact integration of the advection step as explained before. Nev- 
ertheless, one may gain in the number of used velocities that are not constrained any longer to be space-filling ones. In 
contrast to point- wise interpolation schemes, this approach can be applied without compromising exact conservation 
laws or equilibrium properties. Also, due to specific properties of the methodology, the resulting TVLBM can operate 
on adaptive meshes, thereby providing a significant boost of geometrical flexibility especially close to the boundaries 
and in the properties of boundary conditions. 
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II. THE TVLBM MODEL 



Our starting point is the continuum sing le time BGK d HI Hi model written as 

d t Mx,t)+ Cl • V/,(aj,t) = - X - (fi(x,t) - fl eq) (x,t)) 



(1) 



where the left hand side represents the streaming of a probability density function, fi(x,t), to find in the space-time 
location (a?,t) a particle whose velocity v = q is suitably chosen as belonging to a discrete set, thus enforcing the 
desired accuracy order. In terms of the probability density function, we can define macroscopic local variables as the 
density (p), velocity (u) or temperature (T) 



P = ^2fi( x i t ) pu = ^2cifi(x,t) ^pT+^-pu 2 
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with the last two equations that can be combined together to give directly the temperature 



pT=^\ Cl -u\ 2 Mx,t). 



The right hand side of equation (pQ) represents a single time relaxation towards a local Maxwellian equilibrium 

f[ eq \x^t) dependent on (x,t) via the local fields p,u and T. In particular, for the purposes of this paper, the 
following Hermite polynomials representation is adopted 



f ( r\p,u,T) 
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where the shorthand notation of Grad for fully symmetric tensors has been used [Tl|, [35|] . The explicit form of the 
equilibrium distribution, from the standard second order in Hermite polynomials up to the fifth order, is reported 
in appendix A. The presence of a single relaxation time in the evolution equation (pQ) is reproducing only unitary 
Prandtl numbers. This pathology may be removed in different ways [251 l3q |. A simple choice may be considered 
the one proposed in a recent paper by Philippi and coworkers [25], where the right hand side of equation (pQ) is 
supplemented with a local term 



dtfi(x,t)+c r Vfi(x,t) 



1 fl eq \*,t) 
r 9 pT 2 



IL(x,t) : (c/ - u)(ci - u) 



(4) 



The second order tensor is defined in terms of the actual fluctuations with respect to the equilibrium distributions, 
fi — f[ eq ^ , in the following way 



n tj { x ,t) = - fr } (x,t))(4 - Ui )(4 - Uj ). 
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FIG. 1: Left: Diagram for finite volume implementation as reported in the papers of Peng and coworkers [28l l29j]. The points 
P, Pi ,P2,P3,P4,Ps, P6, P7,Ps stand for mesh points while A,B,C,D,E,F,G,H constitute the edges of the control volume where 
the integration of the lattice Boltzmann model is performed. The grey area is where the algorithm is sketched in equations J7J), 
dH} and J9}. Right: A regular (rectangular) mesh with variable spacings Axi (i = l...N x ) and Ayj (j = l...N y ). 



The very general scheme we will consider is therefore 

d t f l (x,t) + c l -Vf l (x,t) = n(x,t) (5) 

with 

Sl(x,t) = -\ (/,(*,*) - f[ eq \x,t)) + l iL^d n(x 9 t) : (Q - u)(q - u). (6) 
III. THE TVLBM EVOLUTION SCHEME 



Following Peng and coworkers [28, 29], the control volume is chosen as a generic polygon ABCDEFGH (see figure 
[1]) surrounding the desired mesh node P so that A, C, E and G are midpoints of edges PPi, PP2, PP3 and PP4 
respectively, while P, D, F and H are the geometric centers of elements PP\P^P2 1 PP2P^P?>iPP?>PiP^ and PPiPgP^ 
respectively. We then treat the polygon ABCDEFGH as made of four elements PABC, PCDF, PEFG, PAHG 
and we focus on the element PABC with the other integrations done in a similar way. The various terms of equation 
(|5|) are then integrated as 



ci • V Jidci =a • tiabIab h Q • n B c^BC — 



PABC z 2 



„ n , (MC) + MP)) , , „ . (8) 

c; • ncpicp 7, 1" c i • n PAiPA ^ 



J PABC 4 

In the above, tiab, Kbc, ncp and rip a are the unit vectors normal to the edges AB, BC CP and PA, and ^s, 
/5c, Zc'P and are the lengths of AP, PC, CP and PA respectively. Finally Spabc represents the surface area of 
the element PABC. 
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It is evident that some of the fluxes over the edges simplify and one can write down explicitly the whole evolution 
involving some weighted combination of the function fi(x,t) in the mesh points surrounding P. For the purposes 
of this paper, we can write down explicitly the evolution for a regular grid of N x x N y rectangular elements with 
edges (Ax)i and (Ay)j. If we consider the physical point x corresponding to the mesh point P = the evolution 

equation over a time lapse dt is 

f l (x,t + dt)=f l (x,t)+dt(n l -c l -Vfi(x,t)), (10) 
and we use the following second order scheme to approximate the various quantities on the right hand side of equation 

m 

c * d * fi « (A*), 1%^ ~ l m + ~ - Te f ^ + h f ^) 

2r y /S S 1 1 1 1 \ 



Mi « (^fii(O) + ^W) + ^(2) + ^fii(3) + ^fij(4) + ^fij(5) + -^(6) + -^(7) + ^^(8)) 

where the notation fi(k) (k = 0, 1, 2, 8) has been used to identify the function t) evaluated in the mesh point 
P/e, which is indeed a first neighbor of P = Pq. 



IV. THE LARGE SCALE LIMIT 



The large scale limit of the previous TVLBM is identified with the compressible Navier- Stokes- Fourier equations 
for an ideal gas and non unitary Prandtl number. The technical procedure leading to such kind of equations is the 
well known Chapman-Enskog expansion [37j which is not detailed in this paper. As a matter of fact, once the kinetic 
equations have been correctly discretized, such a kind of procedure exactly follows the same steps of other calculations 
presented in the literature, both for fully continuum and lattice kinetic equations [H, [HI, l38M40j . The final results 
are summarized in the following set of equations 

dtp + di(fmi) = (11) 

2 

p(d t Ui + UjdjUi) = -diP + dj[r)(diUj + djUi - —Sij(d k u k ))} (12) 

2 

p(d t T + UjdjT) = -Pd k u k + dj(2ndjT) + f](d i u j )(d i u j + djUi - —Sij(d k u k )) (13) 
where we have defined the pressure 

P = pT 

and where the transport coefficients are given by 

7~7~ 77 

rj = c v pT- — — ; v= -\ K = c P pTr 
2r + r g p 

with the specific heats at constant volume and pressure given by 

D D 
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FIG. 2: Left: time history for the transverse velocity in the shear wave decay experiment in a given location y = Right: 
time history for the density in the thermal diffusion mode experiment. In both cases, the results of numerical simulations with 
TVLBM model given in (4|) and ([TO]) are obtained with different relaxation times as reported in the figure. The corresponding 
analytical prediction extracted from the linearized hydrodynamic equations ([11H13P is also reported (solid line) . The details for 
the initial conditions and other simulation parameters are reported in the text. 



V. MEASURING THE TRANSPORT COEFFICIENTS: VISCOSITY AND THERMAL DIFFUSIVITY 

The first numerical experiment we discuss is designed to measure the transport coefficients, i.e. viscosity and 
thermal diffusivity, thus verifying the correct convergence towards the hydrodynamical manifold (jl lti!3j) . To avoid 
the complications of boundary conditions, we choose to study the evolution of the thermal diffusion and transverse 
shear wave modes of linearized hydrodynamics: we expect to see a wave mode decaying as e~ Kk 1 in thermal diffusion 
problems, whereas that of a transverse shear decays as e~ vk £ . All simulations are performed on a L x x L y = 0.4 x 500 
domain, with N x x N y = 2 x 150 grid points and dt = 0.005. Three different choices of the relaxation parameters for 
each transport coefficient experiment are adopted: a) r = 0.2 and r g = 0.1 b) r = 0.2 and r g = 0.2 c) r = 0.2 and 
T g = 0.3 for the transverse shear mode decay and a) r = 0.1 and r g = 0.1 b) r = 0.2 and r g = 0.1 c) r = 0.3 and 
T g = 0.1 for the diffusion mode decay. For the thermal diffusion mode, the initial condition is p = po + e sm(2y7r / Ly) , 
u = 0, and a constant pressure. For the shear wave, the initial condition is u x (x,y) = e sm(2iry / L y ) , u y (x,y) = 0, 
p = po. For all simulated cases, we choose po = I and the perturbation magnitudes are set to e = 0.01 to ensure that 
we stay in the linear regime without influence of the non linear terms. The lattice velocity model used is the D2Q21 
[32| ensuring isotropy up to the eighth order tensors with 21 off-grid velocities. In both cases, the time histories 
of the perturbation magnitudes for velocity (in the kinematic viscosity measurements) and density (in the thermal 
diffusivity measurements) are measured and reported in figure [2j The predicted analytical behaviour is found to be 
well reproduced by the numerical simulations. This is a clear indication that the hydrodynamic equations are very well 
reproduced even with non unitary Prandtl number, as in the numerical simulations we have kept fixed one relaxation 
time of the model and varied the other. 



VI. A TEST FOR COMPRESSIBILITY WITH ENERGY CONSERVATION: THE SHOCK TUBE 

The next numerical experiment to reveal the correct compressible thermohydrodynamical evolution is the one- 
dimensional Sod-Riemann problem [4lj |. We have chosen a two dimensional domain L x x L y = N x dx x N y dy with 
N x = 1500, N y = 2 and dx = dy = 0.0006 so that L x = 1 and L y is so small to make the whole setup result in a 
one dimensional problem. Initially, the gas is at rest (u x = 0) with different states on the two sides of the domain's 
middle point: for x < L x /2 we have set p = pi = 1.0 and P = Pi = 1.0 in LB units while, for x > L x /2, we have set 
P — pr — 0.125 and P = P r = 0.1. This kind of initialization is imposed in the numerics with a very sharp hyperbolic 
tangent profile, ~ tanh((x — L x /2)/£), separating the two half regions x < L x /2 and x > L x /2, with £ = 0.007. In 
figure [3] we compare the numerical solution with the exact theoretical prediction obtained by directly integrating the 
thermohydrodynamical evolution for an inviscid fluid using a finite difference Lax scheme. The solution coming from 
TVLBM is affected by viscous and thermal dissipation but we have chosen a very small relaxation time r = 0.005 so 
as to ensure that, in the observed time lag, the viscous effect has negligible influence. For all the simulation we have 
used a single time relaxation model (j4]) with r g ^> 1 and dt = 10 -6 . Again, the lattice velocity model used is the 
D2Q21 ensuring isotropy up to the eighth order tensors with 21 off-grid velocities. At the edge of the segment L x we 
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FIG. 3: Simulation of the Sod-Riemann shock tube. For two given times (t — 0.1 and t — 0.2) the numerical results obtained 
with TVLBM are compared with the solution of the inviscid Euler equations for a compressible gas (solid line). The numerical 
simulation with TVLBM is done with 21 off- lattice velocities whose properties are reported in [32(]. All the other simulation 
details are reported in the text. 



have set adiabatic boundary conditions (i.e. zero gradient) for all kinetic populations, i.e. fi(—dx,t) = fi(0,t) and 
f l (L x ,t) = fi(L x + dx,t). 

As it is clear from the figure, a rarefaction wave at left, a shock at right and a contact discontinuity at middle are 
observed. In terms of constant regions and positions of discontinuity, the numerical solution agrees very well with the 
theoretical solution for an ideal gas, i.e. with zero viscosity and heat diffusivity. The smooth effect at the contact 
discontinuity and the two ends of the rarefaction wave in the numerical solution results from small but finite viscosity 
and heat diffusivity. When the viscosity and the heat diffusivity are reduced, the contact discontinuity will become 
obviously sharper. 



VII. BOUNDARY CONDITION: DIFFUSE SCATTERING KERNEL 



In the previous sections we tested the algorithm against controlled problems of thermohydrodynamics to benchmark 
the correct convergence towards the Navier-Stokes-Fourier equations for unbounded fluids. The next point in order 
is the investigation of boundary conditions for velocity and temperature fields in the numerical simulations. We will 
first detail the implementation for boundary conditions inspired by the diffuse reflection concept in the rarefied theory 
of gases [H, [39|, |42-45], with benchmarks against known results for the resulting velocity slip and temperature jump. 
Second, we will discuss some original ideas to implement boundary conditions in hydrodynamical problems without 
the emergence of velocity slip or temperature jump at the walls. 

Let us first detail the implementation of the algorithm due to the boundaries. With reference to figure [TJ if we think 
the wall to be located on the node P, the polygons PAHG and PEFG have not to be included in the evolution 
scheme. The flux terms over the edges PA and EP, that are basically omitted in the bulk flow implementation, are 
now taken into account. Given a boundary point, say xw, the resulting scheme is given by 



fi(x w ,t + dt) = fi(x w ,t) +dt(Qi - ci ■ Vfi(x w ,t)) ci • n < 



(14) 
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for all those ingoing populations, i.e. those c\ such that C\ • n < 0. As for the outgoing populations, i.e. those q 
such that C\ • n > 0, we implement a boundary condition inspired by the diffuse reflection concept: the distribution 
functions directed to the walls mix themselves and thermalize to a local Maxwellian before getting reflected into the 
fluid. Before advancing the ingoing populations with (fl4|) . we impose the condition 

fj(xw f \ El, cz . w <0 ^ ' n l M X W,t) f (eq) ( (eg) T (eq) ) 

2^l,c v n>0 \ C l n \Jl \ a w ,-Lw ) 

for the outgoing populations. In this way, the normal velocity to the wall at time t is always zero, and, when the 
system is propagated from t to t + dt, we will have an inward flux of mass which exactly equals the outgoing one, 
i.e. the net gain of mass due to the boundaries is zero. The requirement that condition (fT5j) is exactly satisfied every 
time before the advancing step implies small local depletion/gain of density. Those tiny variations, if necessary, may 
be balanced upon redefinition of the rest population fo . As for the details of the computational scheme, we use the 
following second order scheme to approximate the various quantities on the right hand side of equation (fT4|) 

* 9 °* * (Ax)i + C (Ax)i_i " + - l f ^) 

Wi * (M% {l fi{2) - l m + ^ (5) + ^ (6) - - ^ (3)) ) 

n, « (^fit(o) + ^fiz(i) + ^n«(2) + ^fiz(3) + ^n«(5) + ^n«(6)) • 



VIII. SLIP AND TEMPERATURE JUMP FOR COUETTE FLOWS 



Given the diffuse boundary conditions ([15]) . one may want to investigate the corresponding slip velocity and tem- 
perature jump developing at the walls. For this purpose, we design two distinct experiments to test separately both 
effects. For the slip flow measurements we have chosen an isothermal Couette flow with zero velocity and unitary 
temperature (u^ = 0.0, T^f^ = 1.0) in the lower wall equilibrium and a finite velocity with unitary temperature 
(u^ q) = 0.01, Ti eq) = 1.0) in the upper wall equilibrium. In the case of the thermal jump simulations we have set 

Uw Q ^ = 0.0, Tw Q ^ = 1.005 and Uw^ = 0.0, Tw^ = 0.995 in the lower and upper walls respectively. We then use 
unitary Prandtl numbers (r g ^> 1) with r G [0.0001 : 0.121]. The simulated Couette flow has been confined in a 
two dimensional geometry L x x L y depending on the value of the relaxation parameter r. In particular, we have 
used L x x L y = 0.02^ x 0.035^ with To = 0.001 and the corresponding dt in the simulations has been set equal to 
dt = 0.000001 . The number of grid points has been kept fixed to N x x N y = 2 x 235. Different velocity sets (well 
detailed in recent papers |ll|, [32|) have been used in the numerical simulations, all of them differing in the accuracy 
of the Hermite polynomials of the equilibrium distribution function: a) D2Q9 model with nine space filling velocities 
b) D2Q12 model with 12 off-grid speeds and third order accuracy c) D2Q21 model with 21 off-grid speeds and fourth 
order accuracy d) D2Q2S model with 28 off-grid speeds and fifth order accuracy. The corresponding results for the 
slip length and temperature jump are reported in figure H] and compared with the analytical prediction coming from 
a perturbative analysis of the BGK model [38|, |42|. In particular, the developed velocity slip and temperature jump 
in the aforementioned numerical experiments have been checked against the prediction 



Vslip/ (^j = 1.43684 r T jump / (^j = 1.84074 r (16) 

where (^^j and (^^j are the slope of the velocity and temperature profiles in the 'Navier-Stokes' region [38|, |42| 

away from the boundary layer [56]. It is evident that TVLBM correctly reproduces the desired slip velocity and 
temperature jump, especially in the limit of small r, where we expect the analytical prediction to work well. The 
importance of higher orders in the equilibrium distribution, especially to get the right temperature jump, can be 
appreciated in the right panel of figure HI 
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FIG. 4: Velocity slip and temperature jump from TVLBM model Q and (|10p with single time relaxation r and diffuse boundary 
conditions (|15|) . Two numerical experiments are designed to compute the slip length and temperature jump emerging at the 
walls from the extrapolation of the profiles away from the boundary layers. Details of the numerical simulations are described in 
the text. The numerical results are then compared with the result expected from a perturbative analysis of the BGK equation 
with such diffusive boundary condition [38, 42] and reported in (|16[) . 



IX. BOUNDARY CONDITIONS: AVOIDING TEMPERATURE AND VELOCITY SLIP 

The previous treatment for the boundary conditions is based on the diffuse-reflection idea and, as also demonstrated 
before, is leading to temperature jump and velocity slip at the boundaries. It is anyhow noted that we may want 
to use the local parameters (u w ,T w ) of the local wall equilibrium fi eq \uw q \Tw Q ^) to exactly impose the measured 
velocity and temperature at the wall, i.e. to prevent slip velocity and temperature jump. To do that, we need to 
impose the very same boundary condition as in the previous section with an equilibrium wall velocity and temperature 
chosen as u^ q) + Su^ q) and T& eq) + 8T^ eq) 

fJ Xu/ t ) - T,l, Crn <ofl( X W,t) M)(JM) + g u (eq) m{eq) , S T^) (17) 

Ji(x w ,t)- {eq) {eq) , q) , q) {eq) Ji K u w H-Ot^ a w H-Oi^ )- 

l^l,c r n>0Jl V a w -\-0U w , ± w -T01 w ) 

The variations Suw^ and STw 9>} are computed with an iterative Newton-Raphson procedure in such a way that the 
measured wall velocity and temperature are the desired ones (see Appendix B). It is computationally found that just 
a few iterations (3 or 4) are enough to precisely set the velocity and temperature to the desired values. 
To benchmark the new boundary conditions, the Couette flow between two infinite plates at different temperatures 
and velocities is simulated using the D2Q21 lattice with fourth order accuracy. We have chosen a two dimensional 
domain L x x L y = N x dx x N y dy = 0.2 x 10 with N x = 2 and N y = 120 and the time step has been set equal to 
dt = 0.005. A relevant parameter in this case is given by the Eckert number Ec = U 2 /c v AT , where U is the velocity 
of the upper wall, c v is the constant volume specific-heat and AT is the temperature difference between the walls. 
The velocity is set to zero in the lower wall and, accordingly with the Eckert number, different from zero in the upper 
wall. The simulations were performed using a fixed Eckert number Ec = 2.0 and variable Prandtl number Pr between 
0.3333 and 1.0. In the numerical simulations, to ensure a constant kinematic viscosity and thermal diffusivity, we 
have rescaled the characteristic times with the local pressure p = pT, i.e. r — » ^ and r g —> Then, r is kept fixed to 
r = 0.1 and r g is varied according to the Prandtl number. Overall, as shown in figure EJ the comparison between the 
numerical results and the corresponding analytical estimates for the thermal Couette flows reveals that the TVLBM 
is able to capture correctly the expected behaviour without temperature jump and velocity slip at the boundaries. 



X. GRID REFINEMENT FOR SIMPLE UNIDIRECTIONAL FLOWS 



In this section we explore the possibility to use TVLBM with variable grid mesh in a very simple and controlled 
problem involving thermal hydrodynamics. We choose a thermal Couette flow between two walls at the same temper- 
ature T = 1 with a shear flow imposed by fixing the velocity of the upper/lower wall to U = ±0.1. The computational 
setup is chosen as a two dimensional one with L x x L y = 0.2 x 10.0 where the streamwise length has been set equal 
to L x = N x dx with N x = 2, dx = 0.1 and periodic boundary conditions along it, whereas the vertical length has been 
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FIG. 5: Temperature profile for the thermal Couette flow. We have defined the normalized temperature Tfa t _ c ^ ld ld , and plotted 
it as a function of the normalized distance from the wall x/L y . The Eckert number is kept fixed to Ec — 2.0 while the Prandtl 
number is varied between Pr — 0.3333 and Pr — 1.0. The corresponding analytical profiles are also shown (solid line). All the 
numerical results have been obtained with ([TO)) and Dirichlet boundary conditions imposed as explained in Section IVIIII In 
the inset we report the shear flow in the velocity field normalized with the wall velocity. As we can see, in both temperature 
and velocity profiles, slip is prevented to emerge at the boundaries. 



covered with N y = 60 grid spacings satisfying 



Vj+i ~ Vj 



1, 2, 3 JV W 



(18) 



Vj 



tanh {P4>j ) 
tanh 



X - i 



(19) 



with P > 1 a parameter determining the degree of non uniformity (see also figure [T]) of the mesh, i.e. the larger is /3 
the higher is the non uniformity. For simplicity, we use a unitary Prandtl number obtained with r g ^> 1 and r = 0.01. 
The time step dt has been chosen equal to dt = 0.005 and the D2Q21 model with 21 off-grid speeds and fourth order 
accuracy has been used. As for the non uniform grid, we have chosen /3 = 2.3 with a resulting grid spacing ranging 
from dy = 0.016631 close to the boundaries up to dy = 0.391 in the middle of the channel. Results are reported in 
figure where the refined numerical profile for the temperature (T r ) is compared with the prediction coming from 
stationary hydrodynamics (jl lti!3j) . The temperature profile from this non uniform grid is also compared with the 
temperature profile (T u ) coming from a uniform grid with spacing dy = 0.1666 at fixed N y . To make it visible the 
effect of refinement, one profile has been shifted uniformly with respect to the other with a quantity ST = 0.001. In 
all the numerical simulations, to ensure a constant kinematic viscosity and thermal diffusivity, we have rescaled the 
characteristic times with the local pressure p = pT ', i.e. r — » ^. 



XI. GRID REFINEMENT IN DEVELOPED RB CONVECTION 



In this section we probe the robustness of the algorithm in some non trivial two dimensional setup where thermal 
fluctuations are present, together with non uniform grid spacings. The setup chosen is two dimensional Rayleigh- 
Benard convection [50|452| between two heated walls with different temperatures above the transition point, where 
convective rolls are present and stationary. As a matter of fact, the use of a volumetric formulation may become a 
valuable choice to investigate turbulent convection where we need to well resolve the boundary layer physics. The use 
of an exact stream and collide structure for thermal lattice Boltzmann codes may cause an error source in determining 
the physical properties of the boundary layer due to the presence of spurious, small, departure from the exact linear 
profile in the mean temperature close to the boundary walls [46] . This departure goes together with the existence of 
small spurious transverse velocity for two-three grid layers close to the wall and are due to the existence of discrete 
velocities which connects up to three layers in the lattice inducing non-local boundary conditions effects. Such effects 



10 



1.003 



Refined TVLBM 
theory 








2 



4 



6 



8 



10 



y/L y 



FIG. 6: Temperature profile for the thermal Couette flow with refined grid. The temperature is plotted as a function of the 
normalized distance from the lower wall. The numerical results have been obtained with (|10p and Dirichlet boundary conditions 
for the hydrodynamic fields, imposed as explained in Section IYHI[ a fixed wall velocity U — ±0.1 on both walls (located at 
and y/Ly = 1) and a unitary Prandtl number have been used. The corresponding analytical profile is also shown (solid line), 
as predicted from stationary hydrodynamics. Also, we have used a non uniform grid normal to the walls, whose details are 
reported in (|18H19[) . In the inset, we highlight the effect of grid refinement close to the lower boundary layer, and compare 
it with the corresponding numerical simulation with the same number of grid points arranged in a uniform way (details are 
reported in the text). To make it a clear distinction between the two profiles, we have shifted the profile with uniform grid 
(T u ) with respect to the one with refined grid (T r ) by a constant ST = 0.001. 



can be annoying for the investigation of highly turbulent regimes, where the boundary layer dynamics becomes crucial 
to drive the correct thermal exchange with the bulk [47j. It is numerically observed that this shortcoming can be 
strongly reduced by moving from LBM algorithms using exact streaming to TVLBM based on finite- volume schemes 
as proposed here. 

In what follows, TVLBM numerical simulations are compared against results obtained using finite difference (FD) 
codes for the incompressible case (full details are reported in [H, EH). The computational setup is chosen as a 
two dimensional box L x x L y — 80 x 40 where the streamwise length has been covered with N x = 32 points with 
periodic boundary conditions, while the vertical length has been covered with N y = 64 points and a non uniform grid 
with details reported in ()18til9|) with /3 = 1.5 . Also, the use of a gravitational acceleration g is needed for thermal 
convection. To do that, we implement a general forcing term in the kinetic equations with its exact representation 
(see for example equation (3.15) in [11]). In figures [8] and [3 we make a one-to-one comparison of TVLBM with FD. 
The TVLBM parameters (temperature difference between cold and hot walls, gravity etc..) have been set in such a 
way to not produce strong compressible effects with the same transport coefficients and convection intensity in both 
codes. In particular, the top/bottom wall temperatures have been set equal to X5 = 0.9 and T u = 1.0, with the 
gravitational acceleration equal to g = 0.0001. The transport coefficients correspond to a unitary Prandtl number 
Pr = 1 and Rayleigh number Ra = 8224. The stationary snapshots of the velocity vector field are reported in figure 
[8] where we see a net satisfactory agreement between the two numerical simulations. Further insight is gained by 
checking the details of the thermohydrodynamical profiles for a fixed x as a function of y in figure [71 The stationary 
profiles are very well superposing, as shown for both temperature and velocity field in the streamwise direction. 



We have discussed a volumetric formulation of lattice Boltzmann for compressible fluid flows with active thermal 
fluctuations (TVLBM). The model has been shown to reproduce correctly the large scale behaviour given by the 
Navier- Stokes- Fourier dynamics with and without boundary conditions. The velocity set has been chosen consistently 
with a Gauss-Hermite quadrature and is not necessarily constrained to be a space filling set, thus reducing in number 
the minimal set needed to obtain the correct hydrodynamic behaviour without compromising exact conservation laws 
or equilibrium properties. Also, due to specific properties of the methodology, the resulting method can easily work on 
adaptive meshes, thereby providing a significant boost of geometrical flexibility. For example, it would be extremely 
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FIG. 7: Stationary velocity vector profiles for Rayleigh- Benard thermal convection. The results from TVLBM (left) are 
compared with those from finite difference methods (right). The simulated system has a size L x x L y — 80 x 40, unitary 
Prandtl number Pr — 1, and Rayleigh number Ra = 8424. Other computational details are given in the text. 




y/L y 

FIG. 8: Temperature and streamwise velocity profiles extracted from the stationary snapshots reported in figure [3 We have 
chosen a fixed x — xq — 21.25 and we have plotted the temperature (main figure) and streamwise velocity profile (inset) from 
TVLBM and finite difference (FD) simulations. 



interesting the study of compressible thermal convection at very high Rayleigh numbers [Hoj , especially close to the 
boundaries, where the properties of the thermal boundary layer need to be well resolved to determine the input of heat 
into the system. At the same time, issues related to the generalization of TVLBM to non ideal gases and multiphase 
fluid flows have not been explored systematically in the literature, and interesting lines of research may be envisaged 
[HO, [H, [H for the future. 

Appendix A 

In this appendix we report the details for the equilibrium distribution with successive approximations, from the 
standard second order approximation up to the fifth order one. Given the space dimensionality D, the various terms 
entering the following definition of the equilibrium 

f^\p, u, T) = oji ^2 ^(p, u, T)H\ n) 

n=0 n - 
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are given by the following expressions 

^\p,u,T)n^ = l(uf-u 2 + (T- l)(c? - D)) (20) 



±a$\p,u,T)H\ 3) = j(uf~ 3« 2 + 3(T - l)(c 2 - D - 2)) (21) 



(22) 



±«<V«,t)*<*> = Mf - 6 t 2 + 3 " 4 + ^(tf " * " 2)(«? " - 2 ) " W)- 

^^(cf - 2(D + 2)c 2 + D(D + 2)) 

i4 5) (p,M,T)^ ; (5) = JL («f - 10z* 2 uf + 15u 4 u,) + - 1) (c?u? — (-D + 6)«f - 3cfzV + (3D + 12)u 2 U| ) + 

|(T - l) 2 («,cf - (2D + 8)c?u, + (-D 2 + 6D + 8)u,) 

(23) 

where we have used ui = u - ci, u 2 = u ■ u, cf = ci ■ Ci. 

Appendix B 

In this appendix we detail the technical issues of the boundary condition based on the combination of diffuse- 
reflection scattering kernel and the Newton- Raphson procedure [55| . We start from the kinetic boundary condition 

Y,i, Cvn >o\ c i- n \fi q '(ul q >,Ti } q >) 
For a given time t, let us define two functions 

F 1 (x w ,t,u w ,T w ;u^\T^) =pu w - J^/jcf, 

(25) 



F 2 {x w ,t,u w ,T w -u^ q \T^) =2pT w + pu 2 w -J^/jcf. 



For prescribed u w and T w , we will find Uw q ^ and to satisfy Fi = F 2 = by means of the iterative Newton- 

Raphson procedure. The functions are expanded in Taylor series with respect to generic variations 5uw 9 ^ and 8Tw Q ^ 
in the wall equilibrium velocity and temperature 

/ F 1 (x w , t, u w , T w ; u { : q) + Su^ q) , T ( w eq) + 5T ( w eq) ) \ 
\F 2 (x w ,t,u w ,T w ; vf£ q) + 5u^ q) , T^ eq) + 8Ti eq) ) ) 

/ F 1 (x w ,t,u w ,T w ;v^ q) ,T^ q) ) \ _ Zi, Crn < \ci-n\Mx w ,t) ( A B \ ( 5u { ^ q) \ (26) 







( su ( : q) \ 




S) 


\ 5Ti eq) i 



V F 2 (x w ,t,u w ,T w ;u^ q) ,Ti eq) ) ) £ J>Ci . n>0 l c * " n\f[ eq) (u^ q) ,Ti eq) ) 
+ 0{{5uW) + 0((5TM) 2 ) + 0((5uM)(6TM)), 
where the coefficients are 



l,Ci -Tl>0 



/V i \( d ti eq) W /V xA 

(eq)\ \2^l,c r n>0\ C l ' n \ \ du^ q) ) ) \J-^hc v n>Q C l J I 



M^J £«,«.»>„ l«-»l/« (e,) 
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B = 




£*,c,.n>0 1^ * U \ 



9f 



yl^i,c v n>o c i h J 



l,ci-n>0 



ci 



n\f\ eq) 



C = 



l,ci -n>0 



' d Ae q y 



du 



Ei, ci .n>0 U \ 



Off 



{^2l,c v n>0 C ffl ^) 



(eq) 



,2 , <>/;■■' 



Neglecting the higher order terms and solving linear simultaneous equations with F\ (xw, u w , T w ; t, Uw Q ^ 



5u { ^ q) ,T^ q) + 5T^ q) ) = F 2 (x w , u w ,T w ; t, u^ q) + 5u^ q) ,T^ q) + 5T^ q) ) = 0, we obtain the corrections 
I SuL eq) \ Yi, Cl . n > Wi-n\ti eq \^ ) ,T^ ) ) f D -B\( F^xw^u^T^u^,^) 



\ STi eq) J {AD-Ba)Z l , Cl . n < \c l -n\Mx w ,t)\-C A J \ F 2 (x w ,t,u w ,T w ;u^ q) ,ri eq) ) 
which are added to the solutions 

U: q) ) =U: q A +su&\ 

V / new V / old 

\ / new V / old 

The process is iterated to convergence. 



(27) 



(28) 
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